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Abstract 

A landfill represents a complex and dynamically evolving structure that can be stochastically perturbed by exogenous factors. Both 
thermodynamic (equilibrium) and time varying (non-steady state) properties of a landfill are affected by spatially heterogenous and 
nonlinear subprocesses that combine with constraining initial and boundary conditions arising from the associated surroundings. 
While multiple approaches have been made to model landfill statistics by incorporating spatially dependent parameters on the one 
hand (data based approach) and continuum dynamical mass-balance equations on the other (equation based modelling), practically 
no attempt has been made to amalgamate these two approaches while also incorporating inherent stochastically induced fluctua¬ 
tions affecting the process overall. In this article, we will implement a minimalist scheme of modelling the time evolution of a 
realistic three dimensional landfill through a reaction-diffusion based approach, focusing on the coupled interactions of four key 
variables - solid mass density, hydrolysed mass density, acetogenic mass density and methanogenic mass density, that themselves 
are stochastically affected by fluctuations, coupled with diffusive relaxation of the individual densities, in ambient sun'oundings. 
Our results indicate that close to the linearly stable limit, the large time steady state properties, arising out of a series of complex 
coupled interactions between the stochastically driven variables, are scarcely affected by the biochemical growth-decay statistics. 
Our results clearly show that an equilibrium landfill structure is primarily determined by the solid and hydrolysed mass densities 
only rendering the other variables as statistically “irrelevant” in this (large time) asymptotic limit. The other major implication of 
incorporation of stochasticity in the landfill evolution dynamics is in the hugely reduced production times of the plants that are now 
approximately 20-30 years instead of the previous deterministic model predictions of 50 years and above. The predictions from 
this stochastic model are in conformity with available experimental observations. 
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1 . Introduction 

Municipal waste management (MSW) has traditionally been a supply chain based facility primarily focused on as¬ 
certaining the most cost effective way of disposing household waste, including bio-waste. The concept of a modern 
landfill though stems from the idea of not only cost optimising bulk bio-waste disposal, but also to recycle the bio¬ 
disposables to convert chemical energy to industrially usable electric energy. From a supply chain perspective, this 
constitutes a feedback architecture where the disposable waste produces usable energy that is then fed back to the 
system itself for self-sustenance of the energy production process while simultaneously trafficking the extra energy 
generated for industrial usage [1]. Such operational management of power production from disposable bio-waste 
fundamentally relies on the engineering novelty that could ensure maximum energy production at minimum bio-filler 
consumption while also maximising the profit generated by appropriate disbursement of the energy through the as¬ 
sociated supply chain network [2]. The success of such an “alternative energy” based industry then is inherently 
determined by the accuracy at which the following two factors can be probabilistically evaluated - the start time of the 
production process and the end time line up to which bulk production can be ensured from a plant. 

Collection rates of the output landfill (methane) gas and associated collection efficiency are pivotal in quantifying 
the quality of a production plant and in future planning deliverables based on such production. Results at real methane 
production sites (methanogenic phase) have shown that the production rate and volume could drastically change 
depending on the nature and quality of clay covers, geosynthetic clay liners and geomembrane composite covers with 
the CH 4 [3] emission rates varying from 2.2 to 10,000 mg/m^/d. Aside of the core landfill engineering, alternative 
(methane) production methods in the form of microbial oxidation have been proposed as a cost efficient measure [4]. 
Numerical models, focusing on the methane production rate with respect to the height dependence of landfill sites 
have supported such observations [5? , 7, 8 ] with additional information such as 99% of the methane gas flow at 
the bottom being oxidized across the 0.8 m soil compost column with bulk oxidation of methane occurring within 
the top 0.2 m. The aspects of municipal solid waste management [9, 10] have been areas of recent research interest, 
especially focusing on landfilling impacts [11] and systems analytics [10, 12] based perspectives. Such statistical 
studies have made extensive use of linear programmming algorithms [12, 13] in analyzing multi-phase mixing of 
leachates. In a recent work [14], landfill gas generation data from residual municipal solid waste (RMSW) have 
been utilized to estimate the anaerobic gas generation rate constants. Without having been explicitly mentioned in 
this article [14], the numbers obtained (0.0347-0.0803 y“') seem clearly to indicate the importance of incorporation 
of stochasticity in the landfill gas related mathematical models, albeit the model applied specifically to aerobically 
stabilised MSW. Order of magnitude estimates made in the context of the United Kingdom landfill data also agree 
with such numbers [15]. The methane production time lines of real plants as indicated by numbers in this article 
estimate time periods between 12.5-33 years [14] that are obtained by inverting these rate constants An accurate 
estimation of methane production time lines from landfill sites have evaded estimations from available deterministic 
theoretical models [16, 17, 18, 19, 20 , 21 ] in which these numbers come grossly overestimated by up to 150% further 
confirming the need for improved theoretical models. 

In all these aforementioned studies, explicit incorporation of stochasticity in the otherwise deterministic dynamics 
could prove useful in analyzing the production process as a function of time. This is where stochastic mathematical 
modelling of the degradation rates of landfilled waste and consequent emergence of the hydro-carbon gases (e.g. 

^Production time periods are roughly equal to the inverse of the gas generation rate constants. 
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methane) from the facility assumes importance. While traditional models [16, 17, 18, 19, 20, 21, 22] have been 
successfully able to predict the correct deterministic core of the processes defining a landfill facility, very few of 
these have made any explicit allusion to stochastic modelling. While some attempts [17, 18, 19, 21] have been made 
to analyze the four phased biosolid—>hydrolysed leachate—»acetogenic compounds—»biogas (methane) production 
process, the incorporation of an explicit stochastic uncertainty could emphasise the role of phase heterogeneity in 
the mathematical model. The descriptions in these models clearly indicate the understanding for the need of such 
additional contributions in the model but related attempts were restricted only to the “mean field” probabilistic model. 

The premise of this article is to bridge this information gap between realistic stochastic fluctuations of variables as 
seen in actual landfill sites and assumed deterministic approximation of the same in theoretical modelling as have been 
done. The target is to ensure that not only qualitative facts concerning the landfill dynamics are correctly accounted 
for but also accurate quantitative estimates of decay times of gas production facilities be estimated from the theoretical 
model. This is the objective of this article and would be studied using well established reaction-diffusion formalism 
as detailed below. 

The article is organised as follows. After the derivation and description of the core stochastic model in the fol¬ 
lowing section (Section 2), the temporal dynamics will be analysed in details (Section 3) where the focus will be on 
autocorrelation functions, the squared terms of which will clone the stochastic temporal dynamics. This then will be 
followed by a conclusion section (Section 4) where a summary of the main results will be drawn. 

2. Materials and Methods 

This is a paper on theoretical modelling of a stochastically forced dynamical process. As indicated already, this 
is a non-reductionist scheme geared toward optimised management of resources related to a waste management site, 
primarily focusing on the linear kinetics of the individual variables - solid mass density (ni^^^), hydrolysed mass 
density (u 2 ^^*), acetogenic mass density (« 3 ^“') and methanogenic mass density (« 4 ^"'^) - together with their mutual 
coupled statistics that are often stochastically perturbed by the surrounding environment as well as through spatiotem- 
poral parametric fluctuations. This will be a non-conventional approach that will address the role of complexity in 
the (hierarchical) biochemical pathways’ network. Unlike its predecessors [16, 17, 18, 19, 20] that have generally 
focused on data based phenomenological models, our continuum model will be structured around the well established 
reaction-diffusion scheme that has so often been successfully employed in addressing problems in biology [23, 24], 
fluid mechanics [25] as well as in optimization of resources [7]. An intrinsic strength of our model is its ability to 
identify and analyse the origin of spatial heterogeneity and its time evolution with respect to interacting variables. 
As will be later seen, this is a key issue in determining the eventual decay times of the variables as has been shown 
elsewhere [26]. A reaction-diffusion model with additive noise but with mutually coupled variables is equivalent to 
a multiplicative noise based time evolution that may lead to certain features of nonlinearity starting from a linear 
deterministic structure [26]. One must also note that the choice of the nature of the noise distribution, that is whether 
it is “white” or “colored”, will quantitatively affect the dynamical process. At time scales close to the equilibrium 
time scale, though, this is not expected to affect the result by more than a few percentage at the most. The other tacit 
assumption, and hence a “weakness”, is the assumption of continuum dynamics. Once again, for all known practical 
purposes, this is not expected to have any major effect on the outcome. 
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3. Theory 

In line with the scheme enunciated in earlier works [20, 21], we introduce a generalised hydrolysing stoichiometric 
growth/decay rate as the solution of first order reaction-diffusion kinetics as follows 


R{t) = Af exp(-kf). 


( 1 ) 


where R{t) - reaction rate at time t (kilogram/year) while A and k are respectively the amplitude and decay 
rates expressed in non-dimensional units. The constant a subjectively characterises the specific landfill concerned 
(o' = 1 represents an earlier model [ 21 ] while other limiting values of a correspond to other models detailed before 
[18, 19, 20]). 

5.7. Reaction-Diffusion Model 

The first-order reaction-diffusion kinetics detailed in this section can now be engrained in the core stochastic 
reaction-diffusion model [26, 27] defined in line with the Langevin formalism [28, 29] as follows: 


dt 


vi 


dW'' 

dx^ 


H- J]i(x,t) 


( 2 a) 


dt 


V2 


dx^ 


H- - Aaf exp(-kaf) -H 772(x, f) 


( 2 b) 


dt 


V3 


Qln^ta) 

dx^ 


H- kh'nA^”'^ -H Aaf exp(-kaf) - A„ff exp(-k„,f) -H rjj,{x, t) 


( 2 c) 


diiA'”'^ drnA^"'^ 

-+ Amf exp(-k,„f) H- «h 2 s (0 + 74 (x, t). ( 2 d) 

In the above description, k;, represents the degradation rate of the solid phase while k'h represents the same for the 
biogas (methane) phase. The A’s are the biochemical growth (decay) rate amplitudes and the k’s represent the corre¬ 
sponding decay time scales. As mentioned earlier, the diffusion terms V^n,- reflect the force of homogeneity working 
against the evident heterogeneous phase mixture that is present in a landfill. Although H 2 S is an unavoidable end 
product, since it comes with low relative percentage contribution (~ 0.05%), we would drop this term for subsequent 
calculations at this level of modelling. 

o2 

Here the respective diffusion terms {-^) (f = 1, 2, 3 and 4; the superscripts have been neglected to simplify 
notation, a convention that we will follow throughout this article henceforth) represent the heterogeneous relaxation 
of each variable (ni(x, f) + « 2 (x, f) ni{x,t) + nn{x,t)) defined through diffusion constants v,’s while the 77 ,’s 
characterise uncorrelated Gaussian white noises (an assumption) in three spatial dimensions corresponding to each 
variable n,, as follows: 


< 77 ,(x, t)rjj{x', f) >- 2Di6^(x - x')S(t - t')5ij. (3) 

The above continuum coupled model represented by equations (2a,2b,2c,2d) abide the ensemble averaged mass- 
balance relation: ^ -h ni^A) + >- j-jjg curly bracket “<>” representing the ensemble average 

over all noise realisations. 
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It must be noted that at this minimalist level, we are neglecting effects from noise cross-correlations (e.g. < 
Tji'ni >- 0), an assumption which implies that the origin of stochasticity in solid mass density will not incite an 
identical stochastic response in the hydrolysed, acetogenic or methanogenic mass densities, a reasonable assumption 
at this stage. 

In what follows, we will calculate the ensemble averaged root-mean-square values of the respective autocorrelation 
functions like in the large time (f ^ oo) steady state equilibrium limit. These autocorrelation 

functions represent the experimentally measured stochastic equivalents of their deterministic counterparts (as in [21], 
only for v, = 0 though). 


4. Results 

As detailed in the previous section 3.1), we will start with equations (2a,2b,2c,2d) and then analyse these equations 
in the k — co Fourier space. 

Starting with the the -equation, we get 


dtii dhix 

-iojni (k, cS) — —vk^n\ — i -i- t/i 
miK co) 


III = 


—;'ci) + Vtk^ + ki. ’ 


followed by 


(4) 


Hi 


Ci{t) -< ni(x, t) * ni(x, f -H t) > 


‘I •‘'“I 


dco e 


In the experimentally viable large time equilibrium limit (f 
-autocorrelation function: 


< «T(k, co)ni{-k, -co) > . (5) 

oo), this leads to the following solution for the 


< «i(x, t) * ni(x, t + t) > 


= 4;r“Di 

1 

8_ 


Jo 

4;r^Di 


Vl 

2 sJVlT 


^2g-(viF+i,,)T 

vik^ + kh 

nkh 

exp(-k,,T) H- —-], for (vi < 1). 

^ yy\h 


(6a) 

(6b) 


The above form given in equation (6b) is based on a complex integration within aforementioned limits of the wave 
vector k. For a more accurate expression valid for all k-limits, we will use the form given in equation (6a). This is the 
formula used in the plots shown later. 

Starting from equation (2b), we get 


—icon 2 {k., to) — — V 2 k ff 2 + kf,n 2 — AJiik, co) + 772(k, co). 


(7) 


where h(k., co) - 


Jink, 


:F(1 -H a)e 


,Xl+o.)arcten(^) 




Without much loss of generality we may use the value a - 1 as in [21] to get /i(k, co) - i , where “b” 

alludes to the celebrated Dirac-Delta function as is widely known in the literature [24, 25]. 
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The above prescription leads to 


n2(k, oj) 

that in turn gives 


_ hmiK oj) _ 

{—ioj + V2k^){—ioj + v\k^ + k;,) 


i ^/^AaS(^o + iki) 772(k, w) 

(oj + ikaX-ioj + V2k^) -ico + V2k^ 


( 8 ) 


< n 2 (k, w)*n 2 (-k, —to) >— 


kh < ?7i(k.fa;) * ?7i(-k,-fa;) 

{aP- + V2^fP){<jP + {v\P + kh)^) 


2nAiX'6{oj + iki) * 6{-tJ + ika) /f 2 (k, aS) * / 72 (-k, -to) 
{cP + ka^){tiP + V2^lP) 

(9) 


As previously, for vi > V 2 , the above equations (9) may be complex integrated around all 5 poles in any of the 
halves of the respective Argand diagram to obtain 


< n2(x, f) * U2(x, f + t) > = 


/“'’‘I 

^CXD 

Jo 


du e < n 2 (k, to) * n 2 (-k, -to) > 


InDikh 


1 


V2(Vi^ 2 + (yj + y2)k^kh 


-(e 


-(vik-+k,,)T 


+ e 


-V2k^T- 


) + 


1 


(vi - V 2 )k^kh 


(e 


-Vzk^T _ ^-(vik-+ki,)T 


V2 


( 10 ) 


In line with derivations for the autocorrelation functions corresponding to variables rii and n 2 , it should be noted 
that for the other two variables, our model defines the variable as the independent one while depends on U 4 . 
Starting from equations (2c) and (2d) and following similar algebra as before, we can now evaluate the corresponding 
autocorrelation functions as follows; 


< n4(x, f) * n4(x, f + t) > = 



J" du e 

AtPDa 

r°° 1 

dk- 


Jo 

An'^D^ 

- 

1 

V4 

'Ipvpr 


< «4(k, u) * «4(-k, -u) > 

j^2^-{v4k^+ki,')T 

v^k^ + kh 

Trkl' e“b-V4)t:,/r 

exp(-^//T) + —1-)==“]’ (^"4 < 1). (11) 

^ ■\lv 4 kh' 


In the limit kP 0, the above correlation function takes the limiting value < n 4 (x, t) * U 4 (x, t + t) > = 

20471^- 

^JT ' 

Using the information from equation (11) above, the autocorrelation function defining the acetogenic decay dy¬ 
namics in the limit kP 0 can be obtained as follows 


< «3(x, t) * n3(x, t + t) > - 


/“'’‘I 


= AiPikh 


3 

dk 


r° 

'^£>4 r 

Jo 


< ?f 3 (k, to) * n 3 (-k, -to) > 

]Pg~((yi-V3)k-+kk)T 

[(V4 -H V2)k^ + kP][{V4 - V3)k2 4- kP] ’’’ V3 


{ku 


J:mi„)],(12) 


where k^ax and kmin are respectively defined as the inverse of the smallest and largest length scales in the problem. 
In the context of the landfill model, kmin is the inverse of the landfill diameter while k^ax is the inverse of the landfill 

6 



















Author/Applied Mathematical Modelling 00 (2016) 1-12 


1 


height. So the difference is a small finite number. 

5. Discussions 

The formulas presented above depict the complete dynamical behavior of the spatially heterogeneous landfill 
mechanism that, in turn, is governed by four sub-processes. In order to present the results quantitatively, in the above 
and all future formulations, we will consider same noise strengths, that is D\ — D 2 — - D 4 - Dq without any 

loss of generality. Using vi = 1.0, V 2 = 0.8, kh - 1.0, Di - 0.1, D 2 - 0.05, U/ ^ 0 as the parameter values, 
the solutions of eqn(6b) and eqn(lO) when plotted gives the following time decay profiles. It is to be noted that an 
inherent part of this conclusion relies on the fact that vi V 2 ; in other words, on a heterogeneous spread through 
diffusional relaxation. 

The appended figures compare the root-mean-squared profiles of all four variables after averaging over all stochas¬ 
tic realisations. More specifically, we plot = V< 'h'(x, f) * n,(x, t + t) > (i= 1,2,3,4) against the time difference 

T. As can be seen, this are Gaussian stationary processes [27] implying that in the large time equilibrium limit 
(f ^ 00 ), the respective autocorrelation functions depend only on the time difference between two specific points of 
measurement and ‘ *' '' *' ' 



Figure 1: The solid line represents the decay profile of the root mean square solid mass density autocorrelation 
function (represented by eqn(6b)) with time. 

Fig 1 and Fig 2 show the time decay of solid and hydrolysed waste that are used to generate the acetogenic and 
methanogenic phases shown in Fig 3 and Fig 4. As is to be expected, the solid waste decay rate is steeper than the 
hydrolysed phase indicating that follow-up (methane) production necessitates a slow build-up leading to the target 
deliverables. 

The comparison of all four plots shown in Fig 5 encapsulates the summary of this theoretical model. While the 

decay rates of the individual phases vary, the final decay lines for all phases converge to the golden number of 20 
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Figure 2; The solid line represents the decay profile of the root mean square hydrolysed density autocorrelation 
function (represented by eqn(lO)) with time. 



Figure 3: The solid line represents the decay profile of the root mean square solid acetogenic density autocorrelation 
function (represented by eqn(12)) with time. 
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Figure 4: The solid line represents the decay profile of the root mean square solid methanogenic density autocorrela¬ 
tion 






Figure 5: Comparison of time variations of the mass densities of all four variables. 
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years for all estimated phases. This number is subject to the parameter values used (indicated above). We have 
tested for other realistic parameter values availed from other publications [14] (although this result is for aerobically 
stabilised MSW only, and hence is an order of magnitude comparison only) to confirm that the decay time line ^ 
always conforms to the time window of 15-30 years which is in line with this reference. An understanding of this 
result can be drawn from the fact that our model simultaneously combines stochastic forcing with diffusive relaxation 
of the different densities, together which contribute to this experimentally agreeable timeline. This result is a huge 
improvement on existing deterministic models [21, 22]. 

6. Conclusions 

The implication of this theoretical analysis goes beyond the estimation of accurate gas production decay times and 
favourable comparisons with experimental data. The major benefit of such a study will be in the development of a 
robust business model in which quantitative dependence of such decay rates with varying system parameters will be 
part of the analytical description henceforth. As is not so very difficult to conceive, landfill site structures and engi¬ 
neering depend on the ambience and country specific facilities that may imply wide variation in parameter values. The 
results presented here incorporate all such provisions, including fluctuating parameter values. A pragmatic underpin¬ 
ning with regard to landfill engineering will be the precise quantitative control of parameters and clear ideas about the 
right parametric regime that will ensure gas production at a specified rate. As these rates may vary between different 
sites, as also with the country concerned, such numerical control would ensure easier and more direct improvement 
of existing landfill engineering frameworks. The work will also help municipalities and city councils to make right 
decision for landfill gas mining and implementing a sustainable landfill gas extraction as well energy recovery project. 
We believe that the relevance of this work can be best availed in association with national and supranational waste 
management policy makers. This may be particularly relevant in the context of a circular economy package, the focus 
of the European policy agenda (for instance), that is always more oriented towards landfill diversion and the promotion 
of other disposal options. 

From the technical perspective, this analysis informs us that our stochastic linearly stable model is primarily 
affected by the time dependent hydrolysable decay rate (Rit)) at finite time scales f < 10 years). For larger times, the 
steady state statistics remains unchanged with respect to changes in R{t). As indicated above, the steady methanogenic 
gas production rate is primarily determined by the production rate of hydrolysed mass starting from solid mass with the 
acetogenic density contributing the least in the process. Unlike the previous models [21, 22], our model satisfies the 
mass-balance relation at an ensemble averaged level and not for all deterministic realisations. This ensures dynamical 
equilibrium for all finite times: < ^(ni +n2 + + n^) >- 0. A conclusion that we draw is that of higher decay rate 

for the solid mass density compared to the other three decay processes that is comparable to existing industrial results 
[30]. In arriving at the plots, although specified fixed values of the noise strengths were used, but the linear stability 
ensured that the qualitative deductions obtained from the respective autocorrelation functions remain unaffected by the 
noise strengths. Overall, it is predicted that the decay times of any of the four phases is shorter than the corresponding 
homogeneous deterministic case studied in previous models [18, 19, 20, 21]. All these above facts are directly related 
to numbers that a landfill engineer may make use of in designing the best landfill facility subjective to the given 


^The Gioannis, et al result is used for order of magnitude comparison only, since ours is an anaerobic model as opposed to the aerobically 
stabilised model of Gioannis, et al. 
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conditions. We should like to add that the model can be further extended to address specific structural remits, related 
to the legal paraphernalia of waste disposal in a country specific manner, and adjusting landfill plant lifelines in 
accordance with allowed range of plant structure parameters. The diffusion rates, which are key components in 
defining the landfill lifeline, are subjective of the size of the landfill; such legislations are often guided by European 
rules. Appropriate adjustment of the boundary conditions of our model could straight away address such components. 

Another non-trivial aspect of this analysis is that the introduction of an explicit uncertainty in the dynamical model 
confirms the fact that heterogeneities (vi V 2 9^ 0 and + v\ + 0) play a vital role in the decay rate statistics. Our 
results predict that in absence of heterogeneity, mass distribution will take a far longer time. This was not so very 
obvious in most of the previous models. 

While retaining the core deterministic dynamics, that happens much to be the same in these models [17,18,19,21], 
our novelty lies in the introduction of two key factors that remained dormant in all these models. Firstly, explicit noise 
(stochastic) terms accompany each of our four dynamical equations representing each phase (biogas, hydrolysed 
leachate, acetogenic phase, methanogenous phase). Secondly, in order to incorporate the natural tendency of any 
physical system to neutralise the presence of any heterogeneity, we have incorporated diffusion terms in each of the 
phases that lead to a more generalised multiphase approach where the stochastically forced phases can mix with each 
other through diffusion, n for the stochastic input; rather it is self-consistently derived from the original dynamics 
represents in equation(2a,2b,2c,2d). 

A reassuring quantitative confirmation of the analysis presented here comes from a comparison with real landfill 
descriptions as presented independently by Gioannis, et al [14]. Inverting the decay constants estimated in this anal¬ 
ysis, we reassuringly arrive at the production time line limitations as being between 12.5-33 years that is perfectly 
in harmony with numbers presented in this article (e.g. Fig 5). This quantitative ramification also indicates the need 
for extending the scope of the present linearly stable model into the more realistic nonlinear regime; in other words 
strategy evaluation in planning waste management [6, 31]. This is the next plan along with a multivariate analysis of 
the outputs from the stochastic model based on a Fokker-Planck structure. 
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